clc;
clear;
%% Read data from .dat file
bypass.case.deg = '60'; %the deg of the model
bypass.case.inlet1 = '5'; %the flowrate of inlet1, (0, 3, 5)
bypass.case.inlet2 = '10'; %the flowrate of inlet2, (5, 10, 15, 20, 25)
mu = 3.96e-3;% In pa*s

 if strcmp(bypass.case.deg,'60')
    bypass.DIR.dirbase = 'C:\Users\67072\OneDrive - The Pennsylvania State University\Matlab file\bypass\60deg';
 elseif strcmp(bypass.case.deg,'30')
    bypass.DIR.dirbase = 'C:\Users\67072\OneDrive - The Pennsylvania State University\Matlab file\bypass\30deg';
 elseif strcmp(bypass.case.deg,'90')
    bypass.DIR.dirbase = 'C:\Users\67072\OneDrive - The Pennsylvania State University\Matlab file\bypass\90deg';
 end

file_name = fullfile(bypass.DIR.dirbase, [bypass.case.deg, 'deg',' inlet1 ', bypass.case.inlet1, ' inlet2 ', bypass.case.inlet2, '.dat']);
bypass = importdata(file_name);
%% define the wall coordinates. Varies by cases
x1 = -5.0745;
y1 = 7.83629;
x2 = 0.716965;
y2 = 4.25096;
x3 = -5.97109;
y3 = 4.25094;
x4 = 1.40645;
y4 = -0.64439;
x5 = 1.4754;
y5 = 13.7659;
x6 = 1.4754;
y6 = 4.73358;
x7 = 1.4754;
y7 = -9.26306;
x8 = 4.57811;
y8 = 13.7659;
x9 = 4.57811;
y9 = -9.26306;
walls = [
    x1, y1, x2, y2;% graft wall top
    x3, y3, x4, y4;% graft wall bottom
    x2, y2, x6, y6;% strange wall in corner
    x4, y4, x7, y7;% vessel wall left bottom
    x5, y5, x6, y6;% vessel wall left up
    x8, y8, x9, y9;% vessel wall right
];
%% find wall points
 x=bypass.data(:,1);
 y=bypass.data(:,2);
 u=bypass.data(:,3);
 v=bypass.data(:,4);
 w=bypass.data(:,5); 

% Identify wall points based on the condition
wall_points = ((u ~= 0 & circshift(u, 1) == 0) | ...
               (u ~= 0 & circshift(u, -1) == 0) | ...
               (v ~= 0 & circshift(v, 1) == 0) | ...
               (v ~= 0 & circshift(v, -1) == 0));

num_points = length(wall_points);
% Add wall_points as a new column to the data array
data_with_wall_points = [x, y, u, v, wall_points];

%% find velocity gradient
% Initialize matrices to store gradients
du_dx = zeros(size(u));
du_dy = zeros(size(u));
dv_dx = zeros(size(u));
dv_dy = zeros(size(u));
% Iterate through points and calculate gradients
for i = 2:num_points
    % Calculate the gradients of u with respect to x and y
    du_dx(i-1) = (u(i)-u(i-1))/(x(i)-x(i-1));
    % Calculate the gradients of v with respect to x and y
    dv_dx(i-1) = -((v(i)-v(i-1))/(x(i)-x(i-1)));
    %dv_dy(i) = (v(i)-v(i-1))/(y(i)-y(i-1));
end

%This is because the data were recorded with all x, and then change to the
%next y. The normal iteration only returns 0 y difference.
range = 437; %may vary with angles
for i = 2:num_points
    for j = 1:range %define search range
        if i+j-1<num_points
            if(x(i+j-1)==x(i-1)&&y(i+j-1)~=y(i-1))

                du_dy(i-1) = ((u(i+j-1)-u(i-1))/(y(i+j-1)-y(i-1)));

            end
        end
    end
end

for i = 2:num_points
    for j = 1:range %define search range
        if i+j-1<num_points
            if(x(i+j-1)==x(i-1)&&y(i+j-1)~=y(i-1))

                dv_dy(i-1) = -((v(i+j-1)-v(i-1))/(y(i+j-1)-y(i-1)));

            end
        end
    end
end
%du = [du_dx, du_dy];
%dv = [dv_dx, dv_dy];

%% calculate normal vectors

normal_vector_coordinate = zeros(size(data_with_wall_points, 1), 4);
    for i = 1:num_points
        normal_vector_coordinate(i, 1:2) = data_with_wall_points(i, 1:2);
        if data_with_wall_points(i,5)==1
        x_point = data_with_wall_points(i, 1);
        y_point = data_with_wall_points(i, 2);

        % Initialize distances array
        distances = zeros(1, size(walls, 1));

        % Calculate distances from the current point to all walls
        for j = 1:size(walls, 1)
            x_closest = walls(j, 1);
            y_closest = walls(j, 2);
            x_far = walls(j, 3);
            y_far = walls(j, 4);

            % Calculate distance to the current wall
            distances(j) = point_to_line_distance(x_point, y_point, x_closest, y_closest, x_far, y_far);
        end

        % Find the index of the closest wall
        [~, closest_wall_index] = min(distances);

        % Get the closest wall coordinates
        x_closest = walls(closest_wall_index, 1);
        y_closest = walls(closest_wall_index, 2);
        x_far = walls(closest_wall_index, 3);
        y_far = walls(closest_wall_index, 4);
    
        % Calculate slope of the closer wall
        m1 = (y_far - y_closest) / (x_far - x_closest);

        % Calculate slope of the normal line
        m_normal = -1 / m1;

        x_normal = 1;
        y_normal = m_normal;
        
        % Optional: Normalize the vector if needed
        normal_vector_magnitude = sqrt(x_normal^2 + y_normal^2);
        x_normal = x_normal / normal_vector_magnitude;
        y_normal = y_normal / normal_vector_magnitude;
        %constrain the normal vector direction
        if x_point<=x2 && x_point>=x1 && y_point>=y2 && y_point<=y1
            x_normal = -abs(x_normal);
            y_normal = -abs(y_normal);
        elseif x_point<x4 && x_point>x3 && y_point>y4 && y_point<y3
            x_normal = abs(x_normal);
            y_normal = abs(y_normal);
        elseif x_point<(x5+0.001) && x_point> (x5-0.001) && y_point>=y6 && y_point<=y5 %must have some error 
            x_normal = abs(x_normal);
            y_normal = abs(y_normal);
        elseif x_point<(x8+0.001) && x_point> (x8-0.001) && y_point>=y9 && y_point<=y8
            x_normal = -abs(x_normal);
            y_normal = -abs(y_normal);
        elseif x_point<(x5+0.001) && x_point>(x5-0.001) && y_point>=y7 && y_point<=y4
            x_normal = abs(x_normal);
            y_normal = abs(y_normal);
        elseif x_point<(x6+0.001) && x_point>(x2-0.001) && y_point>=y2 && y_point<=y6
            x_normal = abs(x_normal);
            y_normal = -abs(y_normal);
        else 
            x_normal = 0;
            y_normal = 0;
        end
        normal_vectors= [x_normal, y_normal];
        normal_vector_coordinate(i, 3:4) = normal_vectors;
        % Display or store the results as needed
        %disp(['Normal Vector at Point (', num2str(x_point), ', ', num2str(y_point), '):']);
        %disp(['u_normal = ', num2str(x_normal)]);
        %disp(['v_normal = ', num2str(y_normal)]);
        %disp('---');
        end
    end
%% calculate the normal velocity gradient at each wall point 

x = data_with_wall_points(:, 1);
y = data_with_wall_points(:, 2);
u = data_with_wall_points(:, 3);
v = data_with_wall_points(:, 4);
isWall = data_with_wall_points(:, 5);

x_normal = normal_vector_coordinate(:, 3);
y_normal = normal_vector_coordinate(:, 4);

% Initialize arrays to store velocity gradient components
du_dx = zeros(size(x));
du_dy = zeros(size(y));

dv_dx = zeros(size(x));
dv_dy = zeros(size(y));
gra_temp = 0.0344;
% Loop through wall points
for i = 1:length(x)
    if isWall(i) == 1
        % Use normal vector as directional vectors
        x_normal_i = x_normal(i);
        y_normal_i = y_normal(i);
        
        % Project a line along the normal vector direction
        x_line = x(i) + linspace(0, 10*x_normal_i, 100); % Adjust the length as needed
        y_line = y(i) + linspace(0, 10*y_normal_i, 100); % Adjust the length as needed
        wall_line_x = x2-x1;
        wall_line_y = y2-y1;
        % Find the intersection point with the wall line
        [x_neighbor, y_neighbor] = polyxpoly(x_line, y_line, wall_line_x, wall_line_y);
        
        % Check if intersection points are found
        if ~isempty(x_neighbor) && ~isempty(y_neighbor)
            % Choose the closest intersection point
            [~, min_index] = min((x_neighbor - x(i)).^2 + (y_neighbor - y(i)).^2);
            
            x_neighbor = x_neighbor(min_index);
            y_neighbor = y_neighbor(min_index);
            
            % Calculate finite differences to approximate velocity gradient
            du_dx(i) = (u(x == x_neighbor & y == y_neighbor) - u(i)) / (x_neighbor - x(i));
            du_dy(i) = (u(x == x_neighbor & y == y_neighbor) - u(i)) / (y_neighbor - y(i));

            dv_dx(i) = (v(x == x_neighbor & y == y_neighbor) - v(i)) / (x_neighbor - x(i));
            dv_dy(i) = (v(x == x_neighbor & y == y_neighbor) - v(i)) / (y_neighbor - y(i));
        else
            disp('Intersection not found!');
        end
    end
end

% Assuming dynamic viscosity is stored in a variable mu
mu = 1.0; % Replace this with your actual viscosity value

% Calculate wall shear stress
tau_w = mu * (du_dx + dv_dy);

% Display or use the wall shear stress values as needed
disp(tau_w);
    %% calculate velocity at the wall tangent direction & wall shear stress

% Create an empty array to store the tangent velocity components
%velocity_projection = zeros(size(u));
wall_shear_stress = zeros(size(u));
%wall_shear_stress_x = zeros(size(u));
%wall_shear_stress_y = zeros(size(u));
% Iterate through points in the data array
for i = 1:num_points
    % Check if the point is a wall point
    if wall_points(i)
        % Extract the velocity components at the current point
        velocity_vector = [u(i); v(i)];
        
        % Extract the normal vector at the current point
        normal_vector = [normal_vector_coordinate(i, 3); normal_vector_coordinate(i, 4)];
        normal_vector = transpose(normal_vector);
        % Calculate the tangent vector by rotating the normal vector by 90 degrees
        %tangent_vector = [-normal_vector(2); normal_vector(1)];

        % Calculate the projection of du and dv
        du_prime_dx_prime = dot([du_dx(i),0],normal_vector)/norm(normal_vector);%du'/dx' = ([du/dx,0] . normalvector)/|normalvector|
        dv_prime_dx_prime = dot([dv_dx(i),0],normal_vector)/norm(normal_vector);
        du_prime_dy_prime = dot([0,du_dy(i)],normal_vector)/norm(normal_vector);
        dv_prime_dy_prime = dot([0,dv_dy(i)],normal_vector)/norm(normal_vector);
        wall_shear_stress(i) = mu*((du_prime_dx_prime) + abs(dv_prime_dx_prime) + ...
            (du_prime_dy_prime) + (dv_prime_dy_prime));

        %if wall_shear_stress(i)>0.00019 || wall_shear_stress(i)<-0.0001 %filter some ironic data points
           %wall_shear_stress(i) = 0;
        %end
    end
end
%normal_vector_coordinate(:,5) = wall_shear_stress;
final_data = [x, y, u, v, wall_shear_stress];%final data with x, y, u, v, wss

%% individual plot
%
figure(1);
quiver(x, y, u, v);%plot original vector field
title ('Vector field');

figure(2);
quiver(x, y, du_dx, du_dy, 5);%plot grident of u with respect to normal vector
title('u gradient');

figure(3);
quiver(x, y, dv_dx, dv_dy, 5);%plot grident of v with respect to normal vector
title('v gradient');

figure(4);
quiver(x, y, normal_vector_coordinate(:,3), normal_vector_coordinate(:,4));%plot normal vectors
title('Normal vectors');
%}
%% subplot
%{
figure(1);
subplot(2,2,1);
quiver(x, y, u, v);%plot original vector field
title ('Vector field');

subplot(2,2,2);
quiver(x, y, du_dx, du_dy, 0);%plot grident of u with respect to normal vector
title('u gradient');

subplot(2,2,3);
quiver(x, y, dv_dx, dv_dy, 0);%plot grident of v with respect to normal vector
title('v gradient');

subplot(2,2,4);
quiver(x, y, normal_vector_coordinate(:,3), normal_vector_coordinate(:,4));%plot normal vectors
title('Normal vectors');
%}
decompose_all_velocities(data_with_wall_points, normal_vector_coordinate(:, 3), normal_vector_coordinate(:, 4));

quiver(x, y, u_normal, v_normal,0);
%% Plot the scatter figure
% Extract the x, y, and total wall shear stress values at wall points
%
figure(5);
x_wall = x(wall_points);
y_wall = y(wall_points);
wall_shear_stress_wall = wall_shear_stress(wall_points);

% Normalize the wall shear stress values for color representation
%normalized_wss = (wall_shear_stress_wall - min(wall_shear_stress_wall)) / (max(wall_shear_stress_wall) - min(wall_shear_stress_wall));

% Define a colormap (e.g., jet) for color representation
colormap('jet');
%Plot a 3D scatter plot with color representation
scatter3(x_wall, y_wall, wall_shear_stress_wall, 50,wall_shear_stress_wall, 'filled');
title('Wall shear stress distribusion');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Total Wall Shear Stress');
colorbar; % Add a colorbar to show the color scale
%% user_defined_functions
function distance = point_to_line_distance(x, y, x1, y1, x2, y2)
    % Calculate distance from point (x, y) to the line defined by (x1, y1) and (x2, y2)
    distance = abs((y2 - y1) * x - (x2 - x1) * y + x2 * y1 - y2 * x1) / sqrt((y2 - y1)^2 + (x2 - x1)^2);
end
% decompose the velocity to normal vector
  function [u_parallel_array, v_parallel_array, u_perpendicular_array, v_perpendicular_array,u_prime,v_prime] = decompose_all_velocities(data_with_wall_points, u_normal, v_normal)
    num_points = size(data_with_wall_points, 1);
    u_parallel_array = zeros(num_points, 1);
    v_parallel_array = zeros(num_points, 1);
    u_perpendicular_array = zeros(num_points, 1);
    v_perpendicular_array = zeros(num_points, 1);

    % Calculate the magnitude of the normal vector
   % normal_vector_magnitude = sqrt(u_normal^2 + v_normal^2);

    % Normalize the normal vector
    %u_normal = u_normal / normal_vector_magnitude;
    %v_normal = v_normal / normal_vector_magnitude;

    for i = 1:num_points
       %if data_with_wall_points(i,5) ==1
        u = data_with_wall_points(i, 3);
        v = data_with_wall_points(i, 4);
      
        % Calculate the dot product of velocity and normal vector
            dot_product = u * u_normal + v * v_normal;
    
            % Calculate components parallel and perpendicular to the normal vector
            u_parallel = dot_product * u_normal;
            v_parallel = dot_product * v_normal;
    
            u_perpendicular = u - u_parallel;
            v_perpendicular = v - v_parallel;
    
            % Store the decomposed velocity components in arrays
            u_parallel_array(i) = u_parallel;
            v_parallel_array(i) = v_parallel;
            u_perpendicular_array(i) = u_perpendicular;
            v_perpendicular_array(i) = v_perpendicular;
            u_prime = v_parallel_array+u_parallel_array;
            v_prime = v_perpendicular_array+u_perpendicular_array;
        %end
    end
end
